/*
 * INSANE - Interactive Structural Analysis Environment
 *
 * Copyright (C) 2003-2005
 * Universidade Federal de Minas Gerais
 * Escola de Engenharia
 * Departamento de Engenharia de Estruturas
 * 
 * Author's email :    insane@dees.ufmg.br
 * Author's website :  http://www.dees.ufmg.br/insane
 * 
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 */
package br.ufmg.dees.insane.solution;

import java.io.Serializable;
import java.util.Iterator;
import java.util.Observable;

import br.ufmg.dees.insane.assembler.Assembler;
import br.ufmg.dees.insane.util.IMatrix;
import br.ufmg.dees.insane.util.IVector;
import br.ufmg.dees.insane.util.LinearEquationSystems;
import br.ufmg.dees.insane.util.LoadCase;
import br.ufmg.dees.insane.util.LoadCombination;

/** 
 * A class defining a
 *
 * @author Fuina, Jamile & Pitangueira, Roque
 * @since March 2005
 */
public class StandardNewtonRapson extends Observable implements Step, Cloneable, Serializable{

    private static final long serialVersionUID = 1L;

    protected Assembler assembler; 
    protected LinearEquationSystems sel;
    protected int numMaxIterations;
    protected int currentIteration;
    protected double currentError;
    protected IVector dx;
    protected IVector x;
    protected IVector initialX;
    protected IVector dxp;
    protected IVector dxq;
    protected IVector p;
    protected IVector r;
    private IMatrix cuu;
    protected IVector f;
    protected IVector unbalancedForce;
    protected double loadFactorTotal;
    protected double initialLoadFactor;
    protected int numEq;
    protected double tolerance;
    protected IteractiveStrategy is;
    protected boolean convResult;
    protected byte convType;
    private int label;
    public static final byte FORCE = 1;
    public static final byte DISPLACEMENT = 2;
    public static final byte BOTH = 3;

    public StandardNewtonRapson() {
	this.loadFactorTotal = 0.0;
	this.sel = new LinearEquationSystems();             
	this.convResult = false;
    }    

    public StandardNewtonRapson(int numMaxIterations, double tolerance, byte convType ) {
	this.numMaxIterations = numMaxIterations;
	this.tolerance = tolerance;
	this.convType = convType;
	this.loadFactorTotal = 0.0;
	this.sel = new LinearEquationSystems();				
	this.convResult = false;
    }    

    public IVector getXP(){	
	dxp = sel.solverX(this.cuu,p);	
	return dxp;
    }

    public int getCurrentIteration(){
	return currentIteration;	
    }

    public int getNumMaxIterations() {
	return numMaxIterations;
    }

    public void setNumMaxIterations(int numMI){
	numMaxIterations = numMI;
    }

    public IVector getUnbalancedX(){
	unbalancedForce.sub(r,this.f);	
	dxq = sel.solverX(this.cuu,unbalancedForce);
	return dxq;
    }

    public void assignStepState(double lambda){
	loadFactorTotal += lambda;
	IVector aux1 = new IVector (this.numEq);
	aux1.scale(lambda,p);
	r.add(r, aux1);
	IVector aux2 = new IVector (numEq);
	aux2.scale(lambda,dxp);
	dx.add(aux2, dxq);
	x.add(x, dx);    
	assembler.setXu(x);	
    }

    public boolean getConvergence(){
	return convResult;
    }

    public void setConvergence(){
	Convergence conv = new Convergence(tolerance);
	switch(convType){
	case FORCE:	{
	    convResult = conv.checkConvergence(unbalancedForce,r);
	    break;
	}
	case DISPLACEMENT: {
	    convResult = conv.checkConvergence(dx,x);
	    break;
	}
	case BOTH: {
	    convResult = conv.checkConvergence(unbalancedForce,r) && conv.checkConvergence(dx,x);
	    break;
	}
	}
	this.currentError=conv.getError();
    }

    public void setNumMaxInterations(int nmi){
	numMaxIterations = nmi;
    }

    public void setAssembler(Assembler ass){
	this.assembler = ass;
    }

    public Assembler getAssembler( ){
	return this.assembler;
    }

    public LinearEquationSystems getSel( ){
	return this.sel;
    }

    public double getTolerance() {
	return tolerance;
    }
    public void setTolerance(double tol){
	tolerance = tol;
    }    

    public void setIteractiveStrategy(IteractiveStrategy itestr){
	is = itestr;
    }

    public IteractiveStrategy getIteractiveStrategy( ){
	return this.is;
    }
    public void execute(){
	this.numEq = assembler.getSizeOfXu();
	this.dxq = new IVector(numEq);
	this.dxp = new IVector(numEq);
	this.dx = new IVector(numEq);
	this.r = new IVector(numEq);
	this.p = new IVector(numEq);
	this.f= new IVector(numEq);
	this.unbalancedForce= new IVector(numEq);
	this.x = this.assembler.getXu();
	this.initialX=this.x;

	this.p = this.getProportionalLoad();
	this.r.scale(this.loadFactorTotal,p);
	this.r.add(this.getConstantLoad());

	dxq.zero();			  	
	currentIteration = 1;    			
	for ( int j = 0; j < numMaxIterations; j++ ){

	    this.cuu = assembler.getIncrementalCuu();
	    this.getXP();
	    if ( currentIteration == 1 ){
		this.assignStepState(is.getPredictor());
		this.f=assembler.getFp();
		this.cuu = assembler.getIncrementalCuu();
		this.getUnbalancedX();
	    }
	    else{
		this.f=assembler.getFp();
		this.cuu = assembler.getIncrementalCuu();
		getUnbalancedX();
	    }
	    assignStepState(is.getCorrector());
	    setConvergence();
	    if (!convResult) {
		this.hasChanged();
		this.setChanged();
		this.notifyObservers();
		currentIteration++;
	    }
	    else {
		this.hasChanged();
		this.setChanged();
		this.notifyObservers();
		this.initialLoadFactor=this.loadFactorTotal;
		this.initialX=this.x;
		return;
	    }
	}
	this.loadFactorTotal=this.initialLoadFactor;
	this.assembler.setXu(this.initialX);
	this.f=assembler.getFp();
    }

    /** Returns the Cup*Xp of this Assembler FemModel.<p>
     * The data is stored in the same order of equations.<br>
     *Ex.: CupXp[0] is related to equation -1.
     *@return The Cup*Xp of this Assembler FemModel.
     */
     private IVector getCupXp()
     {
	 IMatrix Cup = assembler.getCup();
	 IVector fed = new IVector(Cup.getNumRow());
	 fed.mul(Cup, assembler.getXp());
	 fed.negate();
	 return (fed);
     };

     /** Returns the PORPORTIONAL Cup*Xp of this Assembler FemModel.<p>
      * The data is stored in the same order of equations.<br>
      *Ex.: CupXp[0] is related to equation -1.
      *@return The Cup*Xp of this Assembler FemModel.
      */
     public IVector getProportionalLoad()
     {	
	 IVector efv = null;
	 IVector aux;
	 LoadCombination lcomb = (LoadCombination)this.assembler.getModel().getLoadCombinationsList().get(0);
	 Iterator ite = lcomb.getLoadCasesList().iterator();
	 assembler.addLoading(lcomb);
	 while (ite.hasNext()){
	     LoadCase lcase = (LoadCase)ite.next();
	     if (lcase.isIncremental()) {
		 efv = this.assembler.getEp();
		 aux = new IVector(this.assembler.getSizeOfXu());
		 aux.sub(this.assembler.getNp(),this.getCupXp());
		 efv.add(aux);
	     }
	 }
	 return (efv);		
     }

     /** Returns the CONSTANT Cup*Xp of this Assembler FemModel.<p>
      * The data is stored in the same order of equations.<br>
      *Ex.: CupXp[0] is related to equation -1.
      *@return The Cup*Xp of this Assembler FemModel.
      */
     public IVector getConstantLoad()
     {
	 IVector efv = new IVector(this.assembler.getSizeOfXu());
	 IVector aux;
	 LoadCombination lcomb = (LoadCombination)this.assembler.getModel().getLoadCombinationsList().get(0);
	 Iterator ite = lcomb.getLoadCasesList().iterator();
	 assembler.addLoading(lcomb);
	 while (ite.hasNext()){
	     LoadCase lcase = (LoadCase)ite.next();
	     if (!lcase.isIncremental()) {
		 efv = this.assembler.getEp();
		 aux = new IVector(this.assembler.getSizeOfXu());
		 aux.sub(this.assembler.getNp(),this.getCupXp());
		 efv.add(aux);
	     }
	 }
	 return (efv);   
     }

     /* (non-Javadoc)
      * @see br.ufmg.dees.insane.solution.Step#getLabel()
      */
     public int getLabel() {
	 return this.label;
     }

     /* (non-Javadoc)
      * @see br.ufmg.dees.insane.solution.Step#setLabel(int)
      */
     public void setLabel(int lb) {
	 this.label = lb;
     }

     /* (non-Javadoc)
      * @see br.ufmg.dees.insane.solution.Step#getLoadFactorTotal()
      */
     public double getLoadFactorTotal() {
	 return loadFactorTotal;
     }

     /* (non-Javadoc)
      * @see br.ufmg.dees.insane.solution.Step#getConvType()
      */
     public byte getConvType() {
	 return convType;
     }

     /* (non-Javadoc)
      * @see br.ufmg.dees.insane.solution.Step#setConvType(byte)
      */
     public void setConvType(byte convType) {
	 this.convType = convType;
     }

     /* (non-Javadoc)
      * @see br.ufmg.dees.insane.solution.Step#getCurrentError()
      */
     public double getCurrentError() {
	 return currentError;
     }

     /* (non-Javadoc)
      * @see br.ufmg.dees.insane.solution.Step#getInitialLoadFactor()
      */
     public double getInitialLoadFactor() {
	 return initialLoadFactor;
     }

     /* (non-Javadoc)
      * @see br.ufmg.dees.insane.solution.Step#setInitialLoadFactor(double)
      */
     public void setInitialLoadFactor(double initialLoadFactor) {
	 this.initialLoadFactor = initialLoadFactor;
     }

     public void setLoadFactorTotal(double loadFactorTotal) {
	 this.loadFactorTotal = loadFactorTotal;
     }
}
